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In  order  to  determine  the  wind  field  in  the  lower  atmosphere  over 

the  entire  globe  one  common  technique  is  to  solve  the  "linear  balance 

equation", 

2      ■*""»■ 

fv  y  +  vy  •  vf  =  g($)  (l) 

where  V   is  the  stress  function  whose  gradient  gives  the  wind  field,  f  is 
the  Coriolis  function 

f  -  2ft  sin  <J>  ,  (2) 

g($)  is  a  known  function  of  the  (known)  geopotential  $,  ft  is  the  earth's 
angular  velocity  and  <j>   is  the  latitude   (polar)  angle  measured  from  the 
equator.   Equation  (1)  is  to  be  solved  on  the  surface  of  a  sphere  of  radius 
a  and  so  is  written 


2 


.      .) 1  a    ,  ,  6  *N  1        6    y|  ,    2ft  cos  4>     d4>  ,..  ,_, 

Sln  * feiT     *J   (cOS   *   6*}   +  — F  -T2 f 2  6*   =   8($)    '         (3) 

(  Y  cos   4>     d  A    )  a 


In  practical  applications  (3)  is  solved  numerically,  by  finite  difference 

methods,  and  the  vanishing  coefficient  of  the  first  term  for  small  <(>  causes 

real  numerical  difficulties.   (The  question  of  whether  (1)  is  an  appropriate 

equation  to  solve  near  the  equator  (for  small  <j>)  is  another  argument  to 

which  we  return  later.)   One  conventional  scheme  to  avoid  these  problems 

is  to  replace  -s — ^   by  a  constant  for  small  <J>  and  in  this  way  continue 

a 
the  solution  of  (3)  for  high  latitudes   to  low  latitudes   by  solving 

2     ~*~         ~* 
f   V  <F  +  Vy  •  Vf  =  g($)  (4) 

for  small  angles  <j>  .   It  is  the  purpose  of  this  note  to  study  the  solutions 
of  homogeneous  versions  of  (1)  and  (4)  in  hopes  of  saying  something  about 
the  solutions  of  the  linear  balance  equation  and  methods  for  solving  it. 


-1- 


We  are  interested  in  the  homogeneous  equation 

sin  <*>  I  -±—    f-  (cos  *  g  )  +  —V  ^4  1  +  cos  *  fl  "  °         (5) 
|  C0S  *  5*        6*      cos2*  6X2  )         5* 

2 
(where  we  have  dropped  the   2fi/a   )  all  the  way  around  the  globe,  so 

we  may  expect  periodicity  on  \    .   We  thus  let 

y  =   eiXn  V(4>)  (6) 

We  also  will  change  the  independent  variable  <j)  by  the  obvious  choice 

x  =  sin<J),  (7) 

and  write  V(<j>)  =  y(x).   Then  (5)  becomes 

«i^  [<>*2>£]  "  ^5  y\    +  (i'-^S  ■  o  (8) 


Guided   by  a  suggestion  in  [1],  we  make  one  further  change  of  variables 

,  n/2 
y(x)  =  (1-x-)    u(x)  (9) 

to  get 

2 

x(l-x2)  ^-§  +  [1  -  (2n+3)x2]  -^  "   (n2+2n)xu  =  0  .  (10) 

,  i  ax 

dx 

Now  (10)  turns  out  to  be  an  equation  studied  almost  100  years  ago  by 

Heun  (see  [1]),  and  is  known  as  Heun's  equation,  or  occasionally  as 

Lame's  equation.   (For  reference  we  list  Babster's  parameters  [1]  for 

(10):   a  =  -1,  b  =  0,  a  =  n,  0  =  n+2,  y  =   1,  6  =  n+1,  e  =  n+1  ).   It  is 

known  that  (10)  has  regular  singular  points  at  x  =  ±  1  and  x  =  °°  ,  and  x  =  0 

Using  the  method  of  Frobenius  [2]  we  can  generate  modified  power  series 

about  each  of  these  points  of  the  form 


OB 

(x-x)    /  ,  c  (x  -  X  ) 
About  the  origin  x  =  0  (the  equator)  the  indicial  equation  for  r  has 
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double  root  of  zero.   About  the  points   x  =  ±  1   (the  North  and  South 
poles)  the  indicial  equation  has  roots   r  =  0,  r  =  -n  . 

The  first  zero  value  for  r  at   x  =  0  means  that  there  is  a  solution 
of  (10)  about  x  =  0  which  is  an  ordinary  power  series.   This  solution 
turns  out  to  be  an  even  function  of   x  and  is  given  by 


u  =  u,(x)  =      Y.         c  x  m  (11) 

1       m^O    m 

=   (m+n/2)(m+l+n/?) 

m+1  /.in2  m 

(m  +  1) 


The  ratio  test  applied  to  (12)  shows  that  (11)  converges  for   |x|  <  1, 
and  the  test  fails  for   |x|  =1.   A  more  subtle  analysis  shows  that  u^  (x) 
converges  for   |x|  =  1.   The  second  zero  value  for  r  at  x  =  0  means 
that  there  is  a  logrtthmic  solution  near  x  =  0,  which  can  be  found  by 
the  method  of  reduction  of  order  [2].   This  solution  is  given  by 

x.   .n+1 
U2(X)=U1(X)  f     (  \  n+lZ   ,    '  <13> 

J      z(l-a)      U^(Z) 

An  analysis  of  (13)  shows  that  u_(x)   does  have  a  logarithmic  singularity 
at  x  =  0,  and,  since  u1 (x)  is  finite  at  x  =  ±  1,   y9(x)   is  singular 
there,  behaving  like   (1  +  x) 

This  is  in  agreement  with  the  indicial  equation  for  (10)  at  x  =  ±  1, 
which  says  that  there  is  one  solution  which  behaves  as  a  constant  (the  r  =  0) 
and  one  which  is  singular  like   (1  ±  x)     .   Thus  the  general  solution 
to  (8)  can  be  written 

y(x)  =  (1  -  x2)n/2  (c1u1(x)  +  c2u2(x))  (14) 
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where  u  (x)   is  finite  everywhere.   Thus  u  (x)   is  singular  at  x  =  0 

and,  at   x  =  ±  1,  is  more  singular  than   (1-x  )    '*•    .   Thus  the  physically 

relevant  version  of   (14)  must  have   c»  =  0,  so  that 

00   inX         n/2 
f(<|»,X)  =  22   e     (cos  <j>)   c  u  (sin  <j>)  (15) 

—  00 

where  the  c  are  the  complex  amplitudes  of  the  longitudinal  harmonics. 

Equation  (15)  is  the  solution  to  the  homogeneous  equation  (5)  and 
is  not  really  useful,  in  itself,  for  solving  the  linear  balance  equation 
(1).   However,  (14)  tells  us  a  great  deat  about  solving  (1)  numerically. 
If  one  started  at  the  poles  or  the  equator   (x  =  ±  1  or  x  =  0)  with  a 
known  (finite)  value  of   4*  then  c        is  necessarily  zero.   Then  taking 
the  exact  solution,   c_   remains  zero.   However,  integrating  numerically, 
one  does  not  have  the  exact  solution,  and  the  roundoff  errors  effectively 
introduce  a  small  amount  of  u_   into  the  solution.   This  would  remain 
small  if  u   were  bounded.   But  integrating  towards  the  equator  (or  the 
poles)  the  unphysical  solution  u~  must  eventually  dominate.   In  other 
words,  no  matter  how  careful  an  integration  scheme  is  used,  because  of  finite 
precision  in  any  computer,  there  is  no  way  to  integrate  the  balance  equa- 
tion (1)  from  the  poles  to  the  equator  and  have  a  realistic  solution  at 
the  equator. 

This  has  been  instinctively  recognized,  a  common  solution  method  for 
the  Northern  Hemisphere  models  {3J  has  been  to  stop  integrating  (1)  at  about 

20  degrees  latitude  and  instead  solve 

2     -*"    + 
2fikV  V  +  Vf  •  Vy  =  g($)  (16) 

o 

where  k  is   sin (20  ).   What  can  one  say  about  the  homogeneous  version 
of  (16)?   Using  the  same  substitutions  as  above  we  obtain,  for  the  homo- 
geneous equation, 


-4- 


(l-x2)un  -  £  [x2  +  2k(n+l)x  -  1]  u'   -  ~(x  +  kn  +  n)u  =  0  (17) 

Equation  (17)  is  not  a  standard  equation,  but  is  amenable  to  a  standard 
power  series  expansion  about  the  origin 

u  -J^x™  (18) 

We  then  obtain  the  recursion  relations 

2a2  =  "  k  al  +  n(n+1>a0 

6a3  =  "  k  a2  +  ("2  +  3n  +  2)  al  +  k  a°  (19) 

m(m-l)a  =  -   (m"1)  a  .   +   n2  +  (2m-3)n  +  (m2-3m-2)  a   „ 
m        k    m-1     L  J  m~2 

(n+m-3) 

k    am-3  ' 

These  involved  equations  show  two  things.   First,  if  k  is  small  (con- 
siderably less  than  1),  they  tend  to  simplify,  and  this  becomes  a  "singular 
perturbation"  difference  equation.  CThis  is  an  almost  unknown  area  which 
the  author  intends  to  investigate  later.)  Secondly,  we  note  that  no  solu- 
tion of  (17)  is  inherently  odd  or  even,  in  contrast  to  both  solutions  of 
(18)  which  are  even  (see  (11)  and  (13)  ).   Thus,  although  any  function  can 
always  be  written  as  the  sum  of  an  odd  and  an  even  function 

f (X)  E  -|  jf (x)  +  f (-X)  J  +  J     jf (X)  -  f (-X)j   , 

numerical  integration  of  (16)  will  always  introduce  some  odd  elements  of 
the  solution.   Since  the  equation  one  is  really  trying  to  solve  is  inherently 
even,  these  portions  of  the  solution  are  unphysical.   Since  all  solutions  of 
(17)  are  bounded  at  x  =  0  >  the  unwanted  portion  of  the  solution  does  not 
grow  drasticly,  so  this  is  not  as  serious  as  trying  to  integrate  (8).   How- 
ever, it  is  expected  that  a  singular  perturbation  analysis  of  the  difference 
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equations  (19)  would  show  that  the  smaller  the  value  of  k  ,  the  more 
closely  the  solution  to  (17)  will  display  the  purely  even  character 
which  (8)  has. 

The  non-homogenous  versions  (1)  and  (16)  could  be  solved  in  power 
series  in  x  =  sin  <f>  and  Fourier  series   in  X.  The  resultant  double  series 
show  exactly  the  features  above.   With  care,  obviously  an  analytic  solution 
in  terms  of  a  double  series  is  possible.   However,  it  certainly  would  be 
difficult  to  get  reasonable  answers  from  such  an  approach.   The  purpose  of 
this  exercise  was  to  demonstrate  that  the  exact  solutions  have  inherent 
problems  and  thus  any  attempt  to  directly  integrate  the  equations  will  give 
difficulty. 

The  author  would  like  to  suggest  an  alternate  approach  to  this  problem 
which  ought  to  be  more  stable,  and  which,  in  other  applications,  is  often 
faster  -  the  finite  element  method.   This  method  has  been  applied  to  some 
elementary  meteorological  problems  [4].   The  basic  idea  is  as  follows: 
instead  of  taking  an  infinite  series  valid  over  the  entire  domain,  one  ex- 
pands  y  as  a  collection  of  low  order  polynomials  -  a  different  polynomial 
for  each  grid  rectangle  in  some  discretization  of  the  globe.   The  coefficients 
for  such  a  ploynomial  are  determined  by  an  integration  of  the  trial  solution, 
times  suitable  functions,  over  the  entire  globe.   The  justification  for  the 

computations  is  either  an  orthogonalization  one  (Galerkin)  or  a  variational 
calculus  one  (Rayleigh-Ritz) .   The  latter  is  more  in  keeping  with  the  use 
of  NVA  (numerical  variational  analysis)  now  used  in  meteorology,  although  it 
will  cause  a  slight  problem  later. 
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For  a  given  function   $  (<f>,A)  it  is  easy  to  see  that  (1)  is  the  Euler 
equation  for 


JJ\f(W)2   +  2g($)Vl 


1(f)-  /  Ijf(VY)   +  2g($)V   a  cos  cf>  d<|)  dA   .  (20) 

D  ' 

where  D  is  the  portion  of  the  globe  of  interest.   If  we  superimpose  on  D 

a  set  of  grid  points  and  connect  these  to  cover  D  with  triangles  (as  opposed 

to  the  rectangles  normally  used  in  finite  difference  solutions  of  (1))  we 

have  a  set  of  domains  D  ,   In  each  of  these  domains  we  express   f  as  a 

i      

simple  function  Y   (a  low  order  polynomial)  whose  coefficients  are  to  be 
i 

determined.   Some  of  the  coefficients  are  determined  by  the  need  for  con- 
tinuity of  Y  and  some  of  its  derivatives  across  the  boundaries  of  D. . 
The  remaining  coefficients  are  determined  by  minimizing  (20) .   The  result 
is  a  matrix  equation  for  these  coefficients  which  is  similar,  but  not 
identical,  to  the  finite  difference  equations  for  (1).   In  particular,  the 
fact  that  one  is  integrating  (20)  over  an  area  means  that  the  vanishing  of 
f   on  a  line  (the  equator)  causes  less  problem  than  before. 

Experiments  run  with  problems  in  stress-strain  mechanics  indicate  that 
using  linear  functions  for  the  Yi,  or  perhaps  quadratics,  will  give  excellent 
accuracy.   The  author  has  not  yet  run  such  an  experiment  for  (20). 

A  question    with  using  (20) (referred  to  before)  is  that  physically 
(1)  may  not  be  the  best  equation, to  use  from  the  poles  to  the  equator. 
Equation  (1)  is  actually 

fV  y     +  Vf  •  VY  =  V2$  (21) 

This  is  really  one  equation  on  2  dependent  variables   ¥  and   $   .   In  the  high 
latitudes   $  is  assumed  known  and  one  solves  for  Y  ,  as  in  (1).   In  the 
equatorial  latitudes  however,  it  is  believed  [5]  that  Y   is  better  known  and 
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one  uses  C21)  to  solve  for  $   .   The.  variational  formulation  of  C21) 
for  $  is 

I  (*)  =  ff         [(V$)2  +  2hCF)$-]a  cos  $d$dA  C22) 

D2 

where  h^)   is  the  left  side  of  (21).   Then  one  could  take  the  domain  for 
(20)  to  be  down  to  some  reasonable  latitude,  say  20  degrees  N,  and  then 

o 

use  (22)  through  the  equator  to   20  S.,  then  use  (20)  down  to  the  South 
pole.   The  difficulty  is  that  V  must  be  precisely  specified  on  the  boundary 
(20°)  for  (20)  and  $   there  likewise  for  (22).   That  is  not  very  realistic. 

Haltiner  has  proposed  [5]  an  iterative  method  to  avoid  the  effects  of 
inaccuracies  in  specifying  the  boundary  conditions.   It  basicly  involves 
solving  three  problems  in  overlapping  regions.   It  involves  solving  for  ¥ 

o  o 

from,  say,  40  N.  to  40  S.,  using  some  reasonable  approximation  (say  $/f) 

for  *F  on  the  boundaries.   This  computation  is  done  by  using  the  observed 

2 
winds  to  find  the  vorticity   £  and  solving  V  ¥  =  £  .   The  resultant  values 

o  o 

of  ¥  at,  say  ,  20  are  used  to  solve  (1)  for  ¥   from  the  poles  to  20 

o         o 

Using  these  values  of   ¥  from,  say,  30   to  20   and  the  values  of  V  from 

O         o  o  o 

±20   to  0   from  step  one,  one  can  then  solve  for   $  from  30  N.  to  30  S. 
using  the  linear  (or  non-linear)  balance  equation.   Because  of  the  inherent 
character  of  the  Laplacean  as  a  smoothing  operator,  this  overlapping  pro- 
cedure should  give  more  reasonable  answers. 

Another  approach  which  minimizes  this  dichotomy  is  to  use  a  finite 
element-Galerkin  approach,  a  non-variational  method.   This  requires  that  we 
use  (21),  the  differential  equation,  rather  than  (20)  or  (22).   One  then 
subdivides  the  entire  globe  into  D.   and  takes   $  and  ¥   to  be  low  order 
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polynomials  in  each  D.   .   In  those  domains  where   $   is  believed  known, 

the  coefficients  of  V.      are  specified.   Then  the  coefficient  of  the  re- 

1 

maining  functions  are  determined  by  insisting  that,  for  all   i   , 
fV  t\>.   +  Viji.  •  Vf  -  V  $ 


If 


D. 

l 


p.  dA  =  0  (23) 


where  the  p.   are  the  low  order  basis  polynomials.   This  is  nearly 
equivalent  to  the  previous  two  integrals,  but  conceptially  places  both 
functions   ¥   and   $  on  the  same  basis.   The  result  of  (23)  is  again  a 
linear  matrix  system  for  the  coefficients  of  the   <J>.   and   ¥ .   ,  the  poly- 
nomial functions  which  represent   $  and  V  in  each  grid  area. 

A  potential  advantage  to  this  approach  is  that  there  is  no  need  what- 
soever to  make  the  grid  spacing  rectangular  or  uniform,  because  there  is  no 

finite  differences  whatsoever.   The   $.   and  ¥ .   are  explicit  functions, 

11 

and  are  differentiated  exactly  by  hand.   By  coverting  to  the  x  ,  X 

coordinate  system  (x  =  sin  <J>) ,  the  coefficients  in  (23)  are  polynomials, 

so  that  the  computations  and  integrations  are  elementary  and  in  fact  are 

exact.   Thus  the  computational  errors  are  only  in  the  use  of  the  $   and 

¥.   ,  and  in  the  data  itself, 
l 
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